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Numerical  Simulation  of  Chemically  Reacting  Flows 
Final  Report  FA9550-12-1-0157 


Overview 

To  an  increasing  extent,  computational  models  are  being  used  in  the  design  of  chemically 
reacting  flow  systems.  This  trend  will  continue  as  models  improve,  computer  power  increases, 
and  the  alternative  of  empirical  testing  becomes  more  expensive.  Of  central  importance  here  is  a 
fundamental  understanding  of  high-temperature,  high-heat-release  systems  —  namely,  flames. 
Such  systems,  often  nonpremixed  and  subject  to  many  competing  physical  effects  involving 
numerous  mixture  components  changing  rapidly  in  space  and  time,  contain  the  essential  features 
inherent  to  practical  engineered  systems.  This  research  looks  to  enhance  the  solution  of  these 
systems  through  a  number  of  computational  methodologies.  In  particular,  we  have  examined 
methods  for  solving  problems  by  constraining  the  temperature  and/or  species  fields;  we  have 
developed  three-dimensional  local  gridding  techniques  and  we  have  utilized  velocity-vorticity 
methods  for  fluid  dynamic  modeling  of  combustion  problems.  Finally,  we  have  combined  high 
order,  high  resolution  spatial  discretization  schemes  with  a  robust  implicit  solution  strategy.  A 
central  premise  of  the  research  discussed  in  this  proposal  is  that  the  advancement  of 
computational  algorithms  and  the  interaction  between  computation  and  experiments  can  play  a 
role  of  equal  or  even  greater  importance,  compared  to  computational  architecture  development, 
in  the  solution  of  these  critical  problems.  All  the  topics  were  designed  to  provide  a  more 
effective  computational/experimental  strategy  for  the  solution  of  problems  of  interest  to  the  Air 
Force. 

Experimentally  Constrained  Computations 

In  our  development  of  computational  models  for  one-dimensional  burner-stabilized  premixed 
laminar  flames,  we  discovered  that  it  was  often  difficult  to  produce  accurate  comparisons  of 
numerical  calculations  with  experimental  data.  This  was  true  even  for  simple  fuels  such  as 
hydrogen  [1].  In  addition,  convergence  difficulties  often  occurred  due  to  the  exponential 
dependence  of  the  temperature  in  the  Arrhenius  chemistry  terms.  As  a  result,  by  specifying  the 
temperature  profile,  we  reduced  the  convergence  difficulties  of  the  Newton-based  solution 
algorithm  and  we  did  not  have  to  model  distributed  heat  losses  (often  radiative).  The  immediate 
benefit  was  that  the  species  comparisons  were  often  dramatically  better  than  when  the 
temperature  was  computed.  This  constrained  approach  to  premixed  flame  computations  is  still 
an  option  in  the  various  PREMIX  versions  that  circulate  in  research  labs  today  [2].  A  variant  of 
the  concept  was  also  utilized  in  a  paper  by  Ashurst  et  ah,  in  which  a  precomputed  velocity 
field  was  utilized  as  an  input  field  in  a  multidimensional  reacting  system  that  approximated  a 
turbulent  reacting  flow  [3].  Generalization  of  these  ideas  has  not  been  explored  in  any  significant 
detail  in  the  solution  of  flames  with  more  complex  chemistry.  As  part  of  our  research  program, 
we  have  utilized  multidimensional  experimental  data  for  the  temperature  and  the  chemical 
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species  to  constrain  the  computation.  This  could  be  in  the  form  of  mole  fractions  or  experimental 
signals.  Several  computational  issues  manifest  themselves  immediately.  Once  the  constrained 
numerical  solution  is  generated,  we  need  to  iterate  the  process  so  as  to  form  revised  experimental 
species  and/or  experimental  signal  information.  As  with  one-dimensional  premixed  flames,  a 
failure  to  converge  implies  an  inconsistency  with  the  model  and  the  experimental  data.  We  have 
applied  this  approach  to  coflow  diffusion  flames  with  simple  fuels  (methane,  ethylene).  A 
comparison  between  experimental  and  computational  soot  volume  fractions  for  a  40%  ethylene- 
air  coflow  diffusion  flame  is  illustrated  in  Figure  1. 


Figure  1.  Soot  profiles  for  a  40%  ethylene-air  coflow  diffusion  flame.  The  contours  on  the  left 
were  generated  using  an  experimentally  measured  temperature  profile.  The  contours  on  the 
right  were  generated  in  a  calculation  where  the  temperature  was  computed. 

Local  Grid  Refinement 

The  local  rectangular  refinement  (LRR)  solution-adaptive  gridding  method,  which  has  been  used 
to  model  numerous  laminar  flames  during  the  past  two  decades,  has  been  extended  to  unsteady 
reacting  applications  in  three  spatial  dimensions  (LRR3D).  Like  LRR,  LRR3D  automatically 
generates  robust  unstructured  grids  by  refining  individual  cells  (not  patches),  uses  multiple-scale 
discretizations,  and  solves  coupled  systems  of  PDEs  via  a  damped  modified  Newton's  method. 
The  grids'  unstructured  nature  produces  a  nonstandard  sparsity  pattern  within  the  Newton's 
method  Jacobian,  as  described  in  the  original  LRR3D  paper  in  2007  [4]. 

Beginning  in  2012,  the  LRR3D  algorithm  was  updated  in  several  ways,  the  main  one  of  which 
was  to  employ  a  full  Jacobian  (not  a  partial  one,  as  in  2007)  and  thus  take  advantage  of  the 
quadratic  convergence  rate  of  Newton's  method.  Furthermore,  the  full  Jacobian  is  implemented 
in  a  memory-efficient  way,  so  that  the  only  Jacobian  blocks  that  are  evaluated  and  stored  are  the 
ones  that  potentially  may  be  nonzero,  based  on  the  local  connectivity  of  points  in  the  grid.  Each 
time  the  grid  adapts,  local  connectivity  is  evaluated  just  once.  This  change  (partial  Jacobian  to 
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full  Jacobian)  was  quite  substantial,  as  it  affected  several  of  the  major  data  structures  in  the  code. 
It  paid  off  in  considerably  faster  convergence  of  the  application  problem  from  the  2007  paper 
(for  one  case,  the  updated  code  required  only  1  Newton  iteration,  down  from  1 1  in  2007). 

More  recent  updates  to  LRR3D  include  the  fact  that  now  many  quantities  appearing  frequently  in 
the  multiple-scale  discretizations  are  precomputed  (each  time  the  grid  adapts)  to  save  CPU  time 
later  during  residual  formation,  and  that  additional  redundancy  has  been  built  into  the  arrays  used 
to  navigate  the  unstructured  LRR3D  grids,  thus  enabling  quicker  access  of  quantities  needed  for 
residual  formation.  These  improvements,  along  with  the  earlier  updates,  have  led  to  a  factor  of 
100  speedup  in  the  grid  adaption  process,  which  enabled  LRR3D  to  be  applied  finally  to 
unsteady  reacting  applications  —  ones  in  which  the  grid  adapts  every  few  time  steps. 
Subsequently,  two  unsteady  reacting  applications  were  investigated  numerically  using  LRR3D, 
as  described  below. 

The  first  application  was  a  3D  extension  of  convection-diffusion-reaction  problem  with  an 
analytical  solution,  whose  purpose  was  to  demonstrate  the  efficiency  and  accuracy  of  LRR3D. 
This  problem  was  governed  by  a  single  nonlinear  PDE.  The  same  problem  was  also  solved 
using  a  structured-grid  code  that  had  discretizations  of  the  same  level  of  accuracy  as  those  in 
LRR3D,  as  well  as  the  same  Newton  solver;  however,  all  data  storage  and  access  in  the 
structured-grid  code  was  greatly  simplified  compared  to  that  in  the  LRR3D  code.  The  structured 
grid  itself  was  equispaced  and  had  the  same  minimum  spacing  as  the  LRR3D  grid.  Conclusions 
were  as  follows: 

•  In  plots  of  spatially  averaged  error  as  a  function  of  time,  LRR3D  results  and  structured- 
grid  results  were  indistinguishable. 

•  LRR3D  used  approximately  530,000  points  (this  value  is  time-averaged,  since  the 
number  of  grid  points  changed  throughout  the  LRR3D  simulation  as  the  grid  adapted), 
while  the  structured-grid  code  employed  approximately  4.2  million  points,  which  was 
fixed  throughout  the  calculation. 

•  LRR3D  required  approximately  one-fourth  of  the  CPU  time  that  the  structured-grid 
calculation  took. 

The  second  application  was  a  3D  extension  of  an  unsteady  2D  solid-solid  alloying  problem 
studied  in  the  mid-1980s  [5].  The  problem  describes  a  reactor  containing  a  mixture  of  solid 
particles  of  aluminum  and  palladium  (Figure  2).  The  particle  radius  varies,  with  larger  particles 
close  to  the  reactor's  central  axis,  smaller  particles  close  to  the  exterior  walls,  and  a  gradual 
variation  in  between.  The  problem's  physics  are  governed  by  two  coupled  nonlinear  PDEs 
involving  two 

unknowns:  temperature  and  solid  fraction.  Values  of  material  properties  appearing  in  the 
equations  are  taken  from  measurements  by  Birnbaum  and  co-workers  [6,7].  Initially  the 
aluminum  and  palladium  are  cold  and  unreacted.  The  bottom  wall  of  the  reactor  is  then 
gradually  heated,  causing  a  diffusion-controlled  alloying  front  to  propagate  through  the  domain. 
The  simulation  terminates  when  the  alloying  is  complete. 

Because  this  application  problem  has  four-fold  symmetry,  only  one-fourth  of  the  physical 
domain  was  modeled.  As  can  be  seen  in  the  figure,  the  alloying  front  (purple)  is  non-planar 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


because  the  local  reaction  rate  depends  on  particle  size.  The  LRR3D  grid  consists  of  the  base 
(equispaced)  grid  plus  four  levels  of  refinement;  the  finest  grid  spacings  automatically  occur  in 
the  region  of  highest  gradients:  surrounding  the  alloying  front.  To  preserve  readability  of  the 
figure,  only  2D  slices  through  the  3D  grid  are  shown.  The  application  was  also  studied  using  the 
same  structured-grid 


Figure  2.  In  the  unsteady  solid-solid  alloying  application,  shown  is  the  computed  solution  at  six  instances 
in  time.  The  alloying  front  is  visualized  by  the  purple  surface,  which  is  where  the  solid  fraction  equals 
0.99  locally.  Isotherms  (other  colors)  are  superimposed  on  four  2D  slices  of  the  LRR3D  grid. 


code  mentioned  above.  For  this  application,  LRR3D's  overall  level  of  accuracy  and  efficiency 
compared  to  the  structured-grid  code  was  similar  to  its  level  of  performance  observed  for  the 
first  application.  Recent  efforts  have  focused  on  deriving  in  3D  Cartesian  coordinates  the 
governing  equations  for  a  vorticity-velocity  formulation  that  includes  non-constant  physical 
properties  (density,  viscosity,  etc.).  Now  that  these  equations  have  been  successfully  derived,  the 
next  challenge  lies  in  developing  multiple-scale  discretizations  of  the  cross-derivative  tenns  — 
in  particular,  discretizations  that  are  as  conservative  as  possible  within  the  constraint  of  using 
only  “allowed”  points  in  the  stencil. 

Vorticity-Velocity  Fluid  Dynamic  Modeling 

Since  2012,  a  novel  vorticity-velocity  formulation  of  the  Navier-Stokes  equations  —  the  Mass- 
Conserving,  Smooth  (MC-Smooth)  vorticity-velocity  formulation  —  has  been  developed.  The 
governing  equations  of  the  MC-Smooth  formulation  include  a  new  second-order  Poisson-like 
elliptic  velocity  equation,  along  with  the  vorticity  transport  equation,  the  energy  conservation 
equation,  and  /VSpcc  species  mass  balance  equations.  Compared  to  the  two  pre-existing  vorticity- 
velocity  fonnulations  (e.g.,  the  Original  (ORIG)  vorticity-velocity  formulation  [8]  and  the 
Modified  (MOD)  vorticity-velocity  formulation  [9]),  the  MC-Smooth  formulation  can  ensure 
mass  conservation  and  solution  smoothness  over  a  broader  range  of  flow  conditions,  and  it 
requires  the  least  CPU  time  to  converge. 

During  the  development  of  the  MC-Smooth  formulation,  we  have  employed  two  axisymmetric 
coflow  laminar  diffusion  flames  as  a  testbed  to  evaluate  the  performance  of  MC-Smooth  and  two 
pre-existing  vorticity-velocity  formulations.  The  first  flame  is  a  confined  flame,  and  the  second 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


flame  is  an  unconfined  flame.  These  two  flames  are  selected  because  they  cover  a  broad 
spectrum  of  operating  conditions  and  the  conclusions  drawn  from  them  are  useful  for  future 
applications. 


Figure  3.  Left:  Normalized  axial  mass  flux  of  the  confined  flame  as  a  function  of  the  axial  position.  At  z  = 
Zmax  =  12.2  cm,  O/O ,  =  0.994  for  the  MC-Smooth  formulation,  Ov/Oi  =  0.984  for  the  modified 
formulation,  and  0/0!  =  0.766  for  the  original  formulation.  Right:  One-dimensional  profiles  of  axial 
velocity  in  the  confined  flame  in  a  portion  of  the  computational  domain  at  z  =  zmax  =  12.2  cm.  At  r  =  2.96 
cm,  vz  =  10.9  cm/s  for  the  MC-Smooth  formulation,  vz  =  1 1.0  cm/s  for  the  modified  formulation,  and  vz  = 
7.0  cm/s  for  the  original  formulation. 

To  quantify  the  mass  conservation  performance  of  each  fonnulation  in  the  confined  flame,  we 
have  computed  the  nonnalized  axial  mass  flux  0/0 1  along  the  axial  direction,  and  the  results  are 
shown  in  the  left  half  of  Figure  3.  As  can  be  seen  from  this  plot,  when  z  =  zmax  =12.2  cm,  the 
MC-Smooth  solution  loses  0.6%  of  its  mass,  the  MOD  solution  loses  1.6%  of  its  mass,  and  the 
ORIG  solution  loses  23.4%  of  its  mass.  This  result  indicates  that  MC-Smooth  and  MOD  can 
ensure  mass  conservation  but  ORIG  cannot.  In  the  right  half  of  Figure  3,  the  corresponding  one¬ 
dimensional  distribution  of  axial  velocity  at  z  =  zmax  =  12.2  cm  are  depicted.  From  this  plot,  it  is 
clearly  observed  that  losing  mass  can  significantly  change  the  prediction  of  the  velocity  field, 
and  that  ensuring  mass  conservation  is  a  crucial  condition  in  obtaining  a  correct  prediction  of  the 
velocity  field. 

For  the  unconfined  flame,  the  two-dimensional  distributions  of  the  axial  velocity  field  predicted 
by  the  MC-Smooth  formulation  and  the  modified  formulation  are  shown  in  the  left  part  of  Figure 
4.  From  this  contour,  we  can  observe  that  the  modified  formulation  is  unable  to  ensure  the 
smoothness  of  the  velocity  field,  and  that  the  nonsmoothness  in  its  velocity  field  will  propagate 
vertically  to  the  downstream.  The  right  part  of  Figure  4  shows  the  one-dimensional  profiles  of 
the  axial  velocity  at  z  =  3.0  cm.  From  this  plot,  it  is  also  observed  that  the  velocity  field  of  the 
modified  formulation  indeed  has  significant  nonsmoothness.  Containing  nonsmoothness  in  the 
velocity  field  is  not  negligible,  because  it  can  further  affect  derived  quantities  such  as  the 
residence  time  and  the  soot  volume  fraction.  The  MC-Smooth  formulation,  on  the  other  hand, 
can  always  ensure  the  smoothness  of  the  velocity  field,  and  this  improvement  is  confirmed  by  the 
two  plots  in  Figure  4.  The  MC-Smooth  fonnulation  also  requires  the  least  CPU  time  to  converge. 
Specifically,  for  the  simulations  related  to  the  confined  flame  (Figure  3),  the  CPU  time  required 
by  MC-Smooth  is  26%  and  16%  lower  than  those  required  by  MOD  and  ORIG,  respectively;  for 
the  simulations  related  to  the  unconfined  flame  (Figure  4),  the  CPU  time  required  by  MC- 
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Smooth  is  39%  and  18%  lower  than  those  required  by  MOD  and  ORIG,  respectively. 
Consequently,  MC-Smooth  not  only  has  a  wider  applicability  than  MOD  and  ORIG,  but  it  also 
requires  a  lower  computational  cost. 


Figure  4.  Left:  Computed  axial  velocity  distribution  of  the  unconfined  flame  in  a  portion  of  the 
computational  domain.  The  left  half  of  the  plot  contains  results  from  the  MC-Smooth  formulation,  and  the 
right  half  of  the  plot  contains  results  from  the  modified  formulation.  Right:  One-dimensional  axial 
velocity  profiles  of  the  unconfined  flame  in  a  portion  of  the  computational  domain.  In  this  plot,  the  solid 
line  is  the  MC-Smooth  prediction,  the  dashed  line  is  the  MOD  prediction,  and  the  dash-dot-dot  line  is  the 
ORIG  prediction.  For  the  data  shown  above,  the  three  lines  mostly  overlap  with  each  other.  The  inset 
highlights  small  differences  among  the  curves. 


In  addition  to  the  improvements  discussed  above,  other  important  advantages  of  the  MC-Smooth 
formulation  include:  (1)  it  does  not  require  the  use  of  a  staggered  grid,  and  (2)  it  does  not  require 
excessive  grid  refinement  to  ensure  mass  conservation.  With  all  these  attractive  features,  the 
MC-Smooth  formulation  is  a  computationally  feasible  approach  that  can  effectively  extend  the 
applicability  of  the  vorticity-velocity  formulation.  Since  2012,  the  MC-Smooth  formulation  has 
been  employed  to  simulate  a  wide  spectrum  of  flames,  including  flames  at  microgravity,  flames 
at  elevated  pressure,  flames  with  soot  formation,  and  flames  with  advanced  radiation  model.  For 
most  flames,  the  MC-Smooth  formulation  can  generate  accurate  predictions  of  the  flame 
structure,  and  very  good  to  excellent  agreement  between  simulations  and  measurements  has  been 
obtained.  In  addition  to  combining  MC-Smooth  with  advanced  numerical  techniques  such  as 
adaptive  grid  refinement,  future  work  will  likely  focus  on  the  application  of  MC-Smooth  to  time- 
dependent  sooting  flames. 

High  Order  Compact  Methods 

Two  key  numerical  challenges  with  compact  methods  have  manifested  themselves  during  recent 
work  on  forced  flames  with  detailed  chemistry,  one  concerning  the  linear  solver,  and  the  other, 
the  Newton  method.  We  have  worked  on  overcoming  these  challenges  by  focusing  on  block- 
based  generalized  saddle  point  preconditioners  and  novel  numerical  techniques  for  compressed 
Jacobian  storage  and  operations.  Ultimately,  a  hybrid  parallelization  strategy  incorporating  both 
task-based  and  domain-based  algorithms  will  make  it  feasible  to  extend  the  validation  of  the 
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solver  to  three-dimensional  problems  while  retaining,  for  the  time  being,  the  global 
approximation  on  which  the  high  order  compact  finite  differences  are  based.  Concurrent  work  in 
the  combustion  diagnostics  portion  of  the  research  has  provided  indispensable  support  for  the 
validation  and  continued  development  of  the  numerical  methods. 

One  critical  issue  that  arose  in  the  C 1  flame  calculations  was  a  form  of  preconditioner  instability 
that  made  it  very  difficult  to  recover  from  a  rejected  time  step  by  simply  reducing  the  stepsize. 
Further  work  using  a  time-dependent  “flame  sheet”  model  considerably  clarified  the  diagnosis  of 
the  problem.  The  relative  simplicity  of  the  flame  sheet  model  has  allowed  us  to  determine  that 
the  preconditioner  instability  is  related  to  certain  features  of  the  fluid  dynamical  subproblem,  in 
particular  (1)  the  presence  of  the  time  derivative  of  density  in  the  continuity  equation,  and  (2)  the 
absence  of  pressure  in  the  same.  Since  the  flames  we  are  studying  bum  in  a  very  low  Mach 
number  regime  where  the  flow  is  hydrodynamically  incompressible,  the  most  natural  set  of 
primitive  variables  includes  the  pressure  at  the  expense  of  the  density,  and  hence  this  time 
derivative  involves  time  derivatives  of  temperature  and  of  all  of  the  mass  fractions  via  the  ideal 
gas  law.  Moreover,  in  a  fully  implicit,  pressure-based  formulation  of  the  problem,  the  pressure 
unknown  is  associated  with  the  continuity  equation,  which  is,  of  course,  independent  of  the 
pressure.  Accordingly,  every  row  in  the  Jacobian  matrix  corresponding  to  the  continuity  equation 
has  not  only  many  strong  off-diagonal  entries  of  size  related  to  the  inverse  of  the  time  step  but 
also  a  zero  on  the  diagonal.  The  effect  of  this  is  potentially  destabilizing  on  any  standard 
incomplete  factorization  algorithm  applied  to  the  low  order  approximate  Jacobian  in  the 
preconditioner  formation  process,  and  it  gets  worse  as  the  time  step  gets  smaller  and  the 
departure  from  diagonal  dominance  increases.  In  some  cases,  such  as  in  our  calculations  of  the 
one-step  flames  with  the  finite-rate  Arrhenius  chemistry  model,  the  preconditioners  constructed 
in  this  way  could  be  made  sufficiently  tractable  for  a  satisfactory  range  of  time  steps  by  various 
stabilization  techniques,  and  the  computations  could  be  completed;  in  others,  it  proved  very 
difficult  or  impossible  to  form  or  maintain  a  sufficiently  stable  and  accurate  preconditioner, 
without  which  the  Newton-Krylov  method  could  not  converge. 

Over  the  past  several  years  we  have  developed  a  detailed  understanding  of  this  difficulty,  as  well 
as  a  plan  for  leveraging  an  important  body  of  recent  (and  ongoing)  research  in  numerical  linear 
algebra  in  order  to  address  it  [10].  The  key  point  is  that  the  discrete  form  of  the  linear  systems 
which  arise  in  the  implicit-compact  solution  process  has  strong  similarities  with  the  form  of  so- 
called  saddle  point  linear  systems,  which  have  received  a  great  deal  of  attention  from  numerical 
analysts  in  the  past  decade  [11].  Typical  examples  of  saddle  point  problems  come  from 
constrained  optimization,  PDE-constrained  optimal  control,  and  incompressible  fluid  dynamics. 
The  common  thread  here  is  the  presence  of  a  constraint.  In  an  incompressible  flow,  the  pressure 
field  is  fixed  by  the  divergence-free  constraint  on  the  velocity  field.  Similarly,  in  the  low  Mach 
number  limit,  the  governing  equations  of  chemically  reacting  flows  describe  a  constrained 
mechanical  system  in  which  the  tiny  hydrodynamic  pressure  can  be  interpreted  as  a  Lagrange 
multiplier  that  imposes  the  proper  divergence  constraint  on  the  velocity  field.  It  can  be  shown 
that  the  only  sources  (S)  of  nonzero  divergence  of  the  velocity  in  an  open  vessel  are  due  to  the 
diffusion  of  heat  and  its  production  by  means  of  the  chemical  reaction,  e.g.,  [12].  Hence,  it  is 
possible  to  formulate  a  low  Mach  number  combustion  problem  as  a  kind  of  incompressible  flow 
problem  with  additional  conservation  equations  for  the  species  and  a  generalized  divergence 
constraint  on  the  velocity  field.  This  insight  is  fundamental  to  some  well-known  splitting 
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methods  for  reacting  flow  problems  [13-15],  yet  the  connection  to  saddle  point  problems  is  not 
clearly  articulated  in  the  combustion  literature,  precisely  because  most  practitioners  of 
computational  combustion  still  avoid  fully  coupled  numerical  methods.  In  the  field  of 
incompressible  fluid  dynamics,  however,  where  implicit  solvers  have  been  used  for  decades,  the 
past  fifteen  years  have  seen  the  advent  of  new,  highly  efficient  saddle  point  preconditioning 
algorithms  for  the  fully  coupled  solution  of  the  steady  Navier-Stokes  equations  [16-19].  Notably, 
these  algorithms  have  also  been  used  in  conjunction  with  a  Newton  solver  [20],  and  recast  in  an 
easily  parallelizable  form  [21].  More  recently,  they  have  been  applied  to  time-dependent 
problems  [22],  and  used  to  study  more  complicated  physics,  such  as  buoyancy  driven  flow  [23]. 
Most  importantly  of  all,  various  generalizations  of  these  algorithms  have  been  devised  for  linear 
systems  which  deviate  from  the  canonical  saddle  point  form  in  one  way  or  another  [24,  25].  The 
time  derivative  of  density  in  the  continuity  equation  is  just  one  element  of  the  reacting  flow 
problem  which  leads  to  problems  with  generalized  saddle  point  structure  (we  note  that  problems 
with  off-diagonal  time  derivatives  have  been  solved  previously  with  numerical  methods  of  the 
kind  we  propose  to  develop  [26]). 

We  have  begun  development  of  these  methods  in  Matlab  [27],  where  the  existence  of  ILUx 
routines  allow  detailed  comparison  with  our  current  algebraic  preconditioning  approach  in  a 
controlled  environment.  The  time-dependent  flame  sheet  problem  provides  a  compelling  test  of 
their  performance  on  a  small-scale  flame  problem.  Applications  to  more  realistic  flames  can 
follow  as  soon  as  suitably  efficient  Schur  complement  approximations  can  be  devised  for 
problems  with  strong,  localized  source  tenns.  One  possibility  is  to  use  block  diagonal 
approximations  for  the  chemistry  terms  (like  in  some  operator  splitting  methods)  and  multigrid 
for  the  flow  variables. 

As  already  mentioned,  the  preliminary  Cl  flame  calculations  undertaken  were  made  possible  by 
replacing  the  JFNK  method  with  a  Newton-like  method  that  employs  a  time-lagged,  low  order 
Jacobian  matrix.  Although  “modified”  Newton  is  an  established  strategy  for  the  efficient 
computation  of  steady  flames  by  “time-dependent”  methods  [28,  29],  our  numerical  experiments 
have  revealed  that  time-lagging  of  the  Jacobian  can  lead  to  inaccurate  results  for  unsteady 
flames,  even  when  the  refresh  rate  of  the  Jacobian  is  relatively  fast  compared  to  the  characteristic 
dynamical  time  scale  of  the  flame  evolution.  Given  that  the  modified  Newton  iteration  is  not  a 
robust  approach  for  true  time-dependent  simulations  and  that  JFNK  has  difficulties  with  small 
chemical  species  which  cannot  be  resolved  well  on  grids  of  modest  size  and  complexity,  better 
Jacobian  options  for  the  nonlinear  solver  will  be  explored  in  the  future.  Automatic 
Differentiation  is  not  a  realistic  option,  since  the  residual  function  that  needs  to  be  differentiated 
is  actually  a  complicated  hierarchy  of  subroutines  that  involves  Chemkin  and  other  third-party 
software  for  the  various  physical  submodels.  A  more  promising  approach  may  be  to  revisit  the 
“Jacobian  component”  methodology,  which  was  first  introduced  in  [30]  and  implemented  in  the 
implicit-compact  primitive  variable  solver  for  nonreacting  flows  that  was  developed  in  the 
earliest  stage  of  our  AFOSR  sponsorship.  This  is  best  understood  as  a  novel  data  compression 
scheme  for  Jacobian  matrices  arising  in  the  nonlinear  solution  of  PDE  problems.  Using  the 
“component  form,”  a  full  (nonsparse)  high  order  Jacobian  can  be  decomposed  into  a  mere  two 
arrays  using  less  memory  than  is  required  to  store  a  nine  blockdiagonal  (sparse)  low  order 
Jacobian.  Moreover,  this  high  order  Jacobian  can  be  applied  to  a  vector  —  the  key  operation  of 
any  iterative  linear  solver  based  on  Krylov  subspaces  —  in  O(N)  operations,  where  N  is  the 
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dimension  of  the  Jacobian.  Although  incredibly  promising  on  paper,  the  effectiveness  of  this 
numerical  approach  was  thought  to  be  dependent  on  the  accurate  analytical  derivation  of  the 
elements  in  the  two  “component”  arrays.  For  nonreacting  flow  problems,  the  robustness  of  this 
process  was  secured  by  an  in-house  code  generator  based  on  the  symbolic  computing  capabilities 
of  Mathematica  [31].  However,  it  proved  extremely  challenging  to  extend  this  software  in  such  a 
way  that  it  could  reliably  derive  the  correct  “Jacobian  components”  for  combustion  problems 
with  detailed  kinetics.  We  believe  that  this  approach  may  be  able  to  be  salvaged,  albeit  at 
marginally  greater  computational  cost  and  with  some  (acceptable)  loss  of  accuracy,  by  replacing 
all  analytically  generated  quantities  with  numerically  estimated  ones.  This  change  should  pose 
no  more  problems  than  the  switch  from  an  analytical  Jacobian  matrix  to  a  numerical  Jacobian 
matrix  —  which,  in  most  normal  circumstances,  is  not  considered  critical  to  the  performance  of  a 
Newton  method  [32]. 
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